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Abstract 

Let f2 be an open, simply connected, and bounded region in R , d > 2, 
and assume its boundary dQ is smooth. Consider solving an elliptic partial 
differential equation -An + yu = f over SI with a Neumann boundary 
condition. The problem is converted to an equivalent elliptic problem over 
the unit ball B, and then a spectral Galerkin method is used to create 
a convergent sequence of multivariate polynomials u n of degree < n that 
is convergent to u. The transformation from Q. to B requires a special 
analytical calculation for its implementation. With sufficiently smooth 
problem parameters, the method is shown to be rapidly convergent. For 
u 6 C 00 (JT) and assuming dfl is a C°° boundary, the convergence of 
\\u — u n \\ Hl to zero is faster than any power of 1/n. Numerical examples 
in R 2 and R 3 show experimentally an exponential rate of convergence. 

1 INTRODUCTION 

Consider solving the Neumann problem for Poisson's equation: 

-Au + 7(s)u = /(s), sen (1) 
du{s) 



g( s ), sedn. (2) 



Assume ft is an open, simply-connected, and bounded region in M. d , d > 2, 
and assume that its boundary 9f2 is several times continuously differentiable. 
Similarly, assume the functions j(s) and f(s) are several times continuously 
differentiable over f2, and assume that g(s) is several times continuously differ- 
entiable over the boundary dfl. 
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There is a rich literature on spectral methods for solving partial differential 
equations. From the more recent literature, we cite [7], [8], [9], and [15]. Their 
bibliographies contain references to earlier papers on spectral methods. The 
present paper is a continuation of the work in [3] in which a spectral method is 
given for a general elliptic equation with a Dirichlet boundary condition. Our 
approach is somewhat different than the standard approaches. We convert the 
partial differential equation to an equivalent problem on the unit disk or unit 
ball, and in the process we are required to work with a more complicated equa- 
tion. Our approach is reminiscent of the use of conformal mappings for planar 
problems. Conformal mappings can be used with our approach when working 
on planar problems, although having a conformal mapping is not necessary. 

In S}2]we assume that ([IJ-© is uniquely solvable, and we present a spectral 
Galerkin method for its solution. In S}3]we extend the method to the problem 
with 7(5) = in f2. The problem is no longer uniquely solvable and we extend 
our spectral method to this case. The implementation of the method is discussed 
in 2] and it is illustrated in <J3J 

2 A spectral method for the uniquely solvable 

case 

We assume the Neumann problem (HJ-© is uniquely solvable. This is true, for 
example, if 

7 (s)>c 7 >0, seH (3) 
for some constant c 7 > 0. For functions u £ H 2 (Q) , v £ H 1 (f2), 



v(s) [-Au(s) + j(s)u] ds = / [ Vit(s) • Vu(s) + j(s)u(s)v(s)] ds 

an on s 



(4) 



Introduce the bilinear functional 

vi,v2) — / [V«i(s) • Vi>2(s) + j{s)vi(s)v2[s)] ds. (5) 
Jn 

The variational form of the Neumann problem (H}-© ^ s as follows: find u such 
that 

A{u,v)=h{v)+i 2 {v), y V £H 1 {Q) (6) 
with the linear functionals defined by 

h(v)= [ v(s)f(s)ds, (7) 

M*)= / v(s)g(s)ds. (8) 
Jan 
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The norms we use for £\ and £2 are the standard operator norms when regarding 
l\ and £2 as linear functionals on H 1 (f2). The functional £\ is bounded easily 
on H 1 (CI), 

\£i(v)\ < ll/MMU* < ll/MMI^- (9) 

Ordinarily, we will use \\v\\i in place of ||w||^fi. 

The functional £2 is bounded (at least for bounded domains fi). To show 
this, begin by noting that the restriction p : — > H 1 / 2 (dfl) is continuous 
[131 Th. 3.37] and the imbedding t : H 1 ' 2 ^®) L 2 (<9ft) is compact [HI Th. 
3.27]. If we further denote by Z g the continuous mapping 

i s : it 1— > / u(s)g(s)ds, u<^L 2 (dVL) 
Jon 

then we see £2 = l g ° t<° p, and therefore £2 is bounded. 
It is straightforward to show A is bounded, 



(10) 



\A(v,w)\ < c A IH^ IH^ , 
c^4 = max {1, hIL}- 
In addition, we assume A is strongly elliptic on H 1 (f2), 

•A («,«)> CeMl VtH 1 ^) (11) 

with some c e > 0. This follows ordinarily from showing the unique solvability of 
the Neumann problem {TJ-C3). If is satisfied, then we can satisfy (fTTj) with 

c e = min {1, c 7 } 

Under our assumptions on A, including the strong ellipticity in (jlip . the Lax- 
Milgram Theorem implies the existence of a unique solution u to ([6]) with 

Mr <i [Pill + INI]. (12) 

Ce 

Our spectral method is defined using polynomial approximations over the 
open unit ball in R d , call it Bd- Introduce a change of variables 

$ : B d i-X n 

onto 

with $ a twice- differ entiable mapping, and let * = : Q ^4 [We 

comment later on the creation of $ for cases in which only the boundary mapping 
4> : dBd — ► c?ri is known.] For v & L 2 (fl), let 

v(x) = v ($ (x)) , xeB d <ZR d 

and conversely, 
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Assuming v £ H 1 (f2), we can show 

V x v(x) = J (xf \7 s v (s) , «=$(z) 
with J (x) the Jacobian matrix for $ over the unit ball Bd, 



J(x) = (£>$) (x) 



<9$ 4 (ir) 



dxi 



X £ B d . 



»J=1 



(13) 



Similarly, 

V s v(s) = K(s) T V x v(x), x = V(s) 
with -ftT(s) the Jacobian matrix for \I/ over SI. Also, 

K(<P(x)) = J(x)~ X . 

Using the change of variables s = $ (sc) , the formula |5]| converts to 

+ ~(($ (x))vi($ (x))v 2 ($ (x)} |det[J(x)]| dx 
{[J(xy T V x vx (x)] T [J (x)~ T V x v 2 (x)} 

I 

+ j(x)vi(x)v2(x)} |det[J(x)]| dx 

{V x vi (%) T A{x)V x V2 (%) + j(x)vi(x)v2(x)} |det [J(x)]| dx 

i 

= A(v 1} v 2 ) (14) 



with 



A(x) ^J(x)~ 1 J(x)~ 



We can also introduce analogues to £\ and l 2 following a change of variables, 
calling them £% and £ 2 and defined on H 1 (Bd). For example, 



v(x)f($(x)) |det[J(x)]| dx. 



We can then convert © to an equivalent problem over H 1 (Bd). The variational 
problem becomes 



A(u,v)=ti(v)+t 2 @), Vv£H 1 (B d ). 



(15) 



The assumptions and results in (j()j)- (fl"Tj) extend to this new problem on H l (Bd). 
The strong ellipticity condition (fTTj) becomes 

A@,v) >c e \\v\\l v£H 1 (B d ), (16) 

|det J(x)\ 



1, max 



xeB d 



\\J(x)\\ 
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where ||J(x)|| 2 denotes the operator matrix 2-norm of J(x) for M. d . Also, 

A(v,w) < caWv^ II x , 

cj[ = < max |det [</(a0]| > max < max || A(x)|| 2 , ||7|| 
Ues d J lx£B d 

For the finite dimensional problem, we want to use the approximating sub- 
space n n = 11^. We want to find u n G H n such that 

A(u n ,v) = hiv) +I 2 (v), We n„. (17) 

The Lax-Milgram Theorem (cf. |H §8.3], [5J §2.7]) implies the existence of u n 
for all n. For the error in this Galerkin method, Cea's Lemma (cf. [H p. 365], 
[3 p. 62]) implies the convergence of u n to u, and moreover, 

< ^ inf ||«-«|k. (18) 

c e »en„ 

It remains to bound the best approximation error on the right side of this 
inequality. 

Ragozin [14] gives bounds on the rate of convergence of best polynomial 
approximation over the unit ball, and these results are extended in [6] to si- 
multaneous approximation of a function and some of its lower order derivatives. 
Assume u S C m+1 (B^)- Using [5J Theorem 1], we have 

inf \u-v i < ^—±w u , m+ i - 19 



with 



w u , m +i (5) = sup sup \D a u(x) - D a u(y)\ 

\ol\— m+1 \|^ — . 



The notation D a u (x) is standard derivative notation with a a multi-integer. 
In particular, for a ~ (ct<i, . . . , ay), 

u ( >~ dx^---dx a / ■ 

When (fT9|) is combined with p8]l . we see that our solutions u„ converge faster 
than any power of 1/n provided u £ C°° {Bdj- 

3 A spectral method for —Aw = / 

Consider the Neumann problem for Poisson's equation: 

-Au=f(s), seQ (20) 
^l=g(s), s edn (21) 



5 



As a reference for this problem, see [5J §5.2]. 

As earlier in ([4]), we have for functions u € H 2 (O) , v € H 1 (f2), 

f f du(s) 

v(s)Au(s) ds = - / Vm(s) ■ Vv(s)ds + / u (s) — — ds (22) 



If this Neumann problem ([2H|) - ([2T|) is solvable, then its solution is not unique: 
any constant added to a solution gives another solution. In addition, if (f20|) - (j2"Tj) 
is solvable, then 



v(s)f(s)ds = / Vu(s) • Vu(s) ds - / v (s) g(s) ds (23) 
Jn Jan 

Choosing v(s) = 1, wc obtain 

f f(s)ds = - f g(s)ds (24) 
Jn JdQ. 

This is a necessary and sufficient condition on the functions / and g in order that 
(|2"0|) - (l2"Tj) be solvable. With this constraint, the Neumann problem is solvable. 
To deal with the non-unique solvability, we look for a solution u satisfying 



u(s)ds = (25) 
Introduce the bilinear functional 

A(vi,v 2 )= / Vui(s) ■ Vv 2 (s)ds (26) 



and the function space 



V = Ui e H 1 (O) : J v(s) ds = o| 



(27) 



A is bounded, 

\A{v,w)\ < \\v\1-l Hiyllj , Wv,w e V. 
From Prop. 5.3.2] A(-, ■) is strongly elliptic on V, satisfying 

A(v,v)>c e \\v\\l, veV 

for some c e > 0. The variational form of the Neumann problem (f20|) - ([2Tj) is as 
follows: find u such that 

A {u, v) = & (v) + l 2 (v) , Vv S V (28) 

with l\ and £2 defined as in (JT])-©. As before, the Lax-Milgram Theorem 
implies the existence of a unique solution u to (|28[) with 

Nli <- Pill + INI]. 
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As in the preceding section, we transform the problem from being defined 
over f2 to being over Bd- Most of the arguments are repeated, and we have 

A{vi,v 2 )= / {VxV^xf A(x)V x v 2 (x)}\det[J(x)}\ dx. 

JB d 

The condition (f2"S")) becomes 

v(x) |det [J(x)]\ dx = 0. 

JB 

We introduce the space 

V = jw e iJ 1 (5) : J u(x)|det[J(a;)]| dx = o|. (29) 

The Neumann problem now has the reformulation 

A (5, v) = It (v) + £ 2 (v) , W e V (30) 
For the finite dimensional approximating problem, we use 

V„ = V n n„ (3i) 

Then we want to find u n £ V n such that 

a (u n ,v)=h (?) + £ 2 (v) , VveV n (32) 

We can invoke the standard results of the Lax-Milgram Theorem and Cea's 
Lemma to obtain the existence of a unique solution u n , and moreover, 

— u n ||i<c inf \\u— v\\\. (33) 

ti£V„ 

for some c > 0. A modification of the argument that led to (fT§|) can be used 
to obtained a similar result for (|3"3"|) . First, however, we discuss the practical 
problem of choosing a basis for V n . 

3.1 Constructing a basis for V n 

Let {(fj : 1 < j < AT^J denote a basis for H n (usually we choose {(fij} to be an 
orthogonal family in the norm of L 2 (Bd))- We assume that <p\(x) is a nonzero 
constant function. Introduce the new basis elements 



C 

with 



tpi =<Pi-^J Vi («) Met [J(x)\ | dx, 1 < j < N% (34) 
C= I |det[J(x)]| dx= \\det[J]\\ Ll (35) 

J B 
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Then tpi = and 

ifj(x)\det[J(x)]\dx= I ipj (x)\det [J (x)]\ dx 

fj(x) |det [J(x)]| dx 



B 



1 

~ C 
= 



|det [J{x)}\ dx 



Thus \(pj : 2 < j < iV„} is a basis of V n and we can use it for our Galerkin 
procedure in (|32|) . 

3.2 The rate of convergence of u n 

Now we estimate inf^y \\u — v\\i\ see (|33|). Recalling ff34|) . we consider the 
linear mapping P : L 2 (Bd) — > L 2 (Bd) given by 

(PS) (x) = u(a;) - i jf |det[ J(y)] | u(y) dy, 
C= ||dct[J]|| L1 ; 
see (j35|) . The mapping P is a projection 

P(PS)(x) = (P2)(a;) - i / |detL%)]| (Pu)(y) dy 



B 

■■u(x)-^J \deb[J(y)]\ u(y) dy- 
I^|det[J(y)]| (s(y) - 1^ |det[J(z)]| u(*)) 

- ^ jf |det[J(y)] | %) d y-^J B |det[J(y)] \ u(y) dy 
i f |det[J(y)]| dy / |det[J(z)]| u{z) dz 



u(x) - ^ J^\det[J(y)}\ u(y) dy 
(Pu)(x) 
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So P 2 = P and P is a projection with ||P||i2_ >i 2 > 1 and 
\\Pu\\ 2 = \\u-± f |det[J(y)]| u(y) dy\\ L , 



< \\u\\ L 2 + 



< \\u\ 



L- 



c 



|det[J(y)]| u{y)dy 
Aet[J]\\ L 2\\u\\ L 2 



111 



L- 



.d/2 



I n*/* ||det[J]|| £2 
r(l + ±d) \\det[J]\\ L1 



r(i + irf) 



( C auchy- Schwarz ) 



CP\\U\\ L 2 



which shows ||P|| L 2^ L 2 < c P and V := P(iP(PM)), see ([29]). For w G iP^Pd) 
we also have Pm G H 1 (Bd) and here we again estimate the norm of P: 



\\Pu\\ 2 m 



\Pu\ 



12 

< 41N2 + IIVulll 



|V(Pu) 



since V(Pu) — Vit. Furthermore cp > 1, so 

HPulliri <4Nli» + 4llv«iii a 

= c 2 P (\\u\\h + \\Vu\\h) 
= c%\\uf H1 
\\Pu\\ H i <c P \\u\\ H i 

and we have also ||P||#i_,#i < cp. For u G V = P(H 1 (B)) we can now estimate 
the minimal approximation error 



mm \\u-p\\ Hl 



_ min \\Pu — p\\hi 
?ev„ 

= min || Pu — PpWfft 
pen,, 

< mm \\P\\ H i^ H i\\u- p\\ Hl 
pen,, 

< cp min ||u — pllijj 

pen„ 



P is a projection 
and u G image(P) 

because V„ = P(n„) 



and now we can apply the results from [5J. 



4 Implementation 

Consider the implementation of the Galerkin method of Sj2] for the Neumann 
problem CO)- ([2]) over Vt, by means of the reformulation in (p~5|) over the unit ball 
Bd- We are to find the function u n G II n satisfying (fTS)) . To do so, we begin by 
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selecting a basis for II„, denoting it by {(fi, . . . , pn}, with N = N n = dimll„. 
Generally we use a basis that is orthonormal in the norm of L 2 (82). It would 
be better probably to use a basis that is orthonormal in the norm of H 1 (Bd); 
for example, see [IB]- We seek 



N 



(36) 



fc=i 



Then (|TT)) is equivalent to 



N 

fe=i 



dip k (x) dip e (x) 

~dx- dx \-i{x)(pk{x)tpe{x) 



|det [J(x)]| dx 



= I f (x) tpe (x) |det [J(x)}\ dx 

JB d 



(37) 



9B d 



g (x) ipe (x) \Jbd y (x)\ dx, 



£=1,...,N 



The function | Jbdy{x)\ arises from the transformation of an integral over <9f2 to 
one over dBd, associated with the change from £2 to £2 as discussed preceding 
(|15p . For example, in one variable the boundary dfl is often represented as a 
mapping 

x(0) = (xi(0), X 2(B)), 0<6<2tt. 
In that case, \Jbdy{x) \ is simply \\' (0)\ and the associated integral is 



g(x(O))w(x(0))\x'(o)\ de 



In p7p we need to calculate the orthonormal polynomials and their first partial 
derivatives; and we also need to approximate the integrals in the linear system. 
For an introduction to the topic of multivariate orthogonal polynomials, see 
Dunkl and Xu [10] and Xu [17]. For multivariate quadrature over the unit ball 
in R d , see Stroud [IB] . 

For the Neumann problem (|20 | - ([2T]) of <J2 the implementation is basically 
the same. The basis {tpi, ■ ■ ■ , ¥>n} is modified as in (|34|) . with the constant C 
of (|35[) approximated using the quadrature in (|42p . given below. 



4.1 The planar case 

The dimension of IL, is 



1 



N n = - (n + 1) (n + 2) 



(38) 



For notation, we replace x with (x, y). How do we choose the orthonormal basis 
{tptix, y)}^! for II n ? Unlike the situation for the single variable case, there are 
many possible orthonormal bases over Bd = D, the unit disk in R 2 . We have 
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chosen one that is particularly convenient for our computations. These are the 
"ridge polynomials" introduced by Logan and Shepp [12] for solving an image 
reconstruction problem. We summarize here the results needed for our work. 
Let 

v„ = {p £ n„ : (P, Q) = VQ G n„_!} 

the polynomials of degree n that are orthogonal to all elements of II„_i. Then 
the dimension of V„ is n + 1; moreover, 

n„ = v e Vi © • • • © v„ (39) 

It is standard to construct orthonormal bases of each V n and to then combine 
them to form an orthonormal basis of II„ using the latter decomposition. As 
an orthonormal basis of V„ we use 

1 7T 

fn,k(x,y) = —;=U n {x cos (kh) +ys'm(kh)) , (x,y) £ D, h= — — (40) 
y/ir n + 1 

for k = 0, 1, . . . , n. The function U n is the Chebyshev polynomial of the second 
kind of degree n: 

UJt) = sm ( n + 1 ) 6 ' ; t = cos6, -l<t<l, n = 0,l,... (41) 
sin 

The family {<^n.fc}^ =0 is an orthonormal basis of V n . As a basis of n„, we order 
{fn,k} lexicographically based on the ordering in (|4TJ)) and (f3"9")) : 

{<Mfcl = {V?0,0, <Pl,0, ^14: ^2,0, ■ ■ ■ , Vn,0, ■ ■ • >^n,n} 

To calculate the first order partial derivatives of p n ,k(%> y), we need U n (t). The 
values of U n (t) and U n (t) are evaluated using the standard triple recursion 
relations 

U n+1 {t)=2tU n {t)-U n - 1 {t) 

U' n+1 (t) = 2U n (t) + 2tU n (t) - C^_ x (t) 

For the numerical approximation of the integrals in (|37[) , which are over B 
being the unit disk, we use the formula 

J g g{x,y) dsdy « £ (r,, Wl -^_ ri (42) 

Here the numbers are the weights of the (g + l)-point Gauss-Legendre quadra- 
ture formula on [0, 1]. Note that 

1 <? 

p(x)dx = y]p(n)uji, 

1=0 

for all single- variable polynomials p(x) with deg (p) < 2q + 1. The formula (|42|) 
uses the trapezoidal rule with 2q + 1 subdivisions for the integration over Bd in 
the azimuthal variable. This quadrature is exact for all polynomials g € Tl2q- 
This formula is also the basis of the hyperinterpolation formula discussed in 

ra- 



il 



4.2 The three dimensional case 



In the three dimensional case the dimension of H n is given by 
and we choose the following orthogonal polynomials on the unit ball 

/ \ (0,m— 27 + 4) no -, \ ri f \ 

y m ,],l3{ x ) = c mjPj m x \\ - l)£>0,m-2j [X) 

= c m A\x\\ m - 2l pT m - 2jH \2\\xf - l)S p , m - V (JL) , (43) 
j = 0,...,[m/2], /3 = 0,l,...,2(m-2j), m = 0,l,...,n 

The constants c mj are given by c m> j = 24 + ~~ : '; and the functions p^'" 1 2j+2 ' ) 
are the normalized Jacobi polynomials. The functions Sp. m -2j are spherical 
harmonic functions and they are orthonormal on the sphere S 2 C M 3 . See 
[TOl [3] for the definition of these functions. In [3] one also finds the quadrature 
methods which we use to approximate the integrals over -Bi(O) in (|14| and (fT"5|) . 
The functional £2 in (jT5]) is given by 

Mv)= r r g ^(T(i,0,4>))) (44) 

J Jo 

■ ||($oT) e (l,M) x (*oTV(l,M)||«(*(T(l,M)))d0d0 

where 

T(p, 0, </>) := p (sin(0) cosf», sin(0) sin (</>), cos(6>)) (45) 

is the usual transformation between spherical and Cartesian coordinates and 
the indices denote the partial derivatives. For the numerical approximation of 
the integral in (|44[) we use traezoidal rules in the <f> direction and Gaufi-Legendre 
formulas for the 9 direction. 



5 Numerical examples 

The construction of our examples is very similar to that given in [3] for the 
Dirichlet problem. Our first two transformations <& have been so chosen that 
we can invert explicitly the mapping to be able to better construct our test 
examples. This is not needed when applying the method; but it simplifies 
the construction of our test cases. Given $, we need to calculate analytically 
the matrix 

A(x) = J(x)" 1 J( X y T . (46) 
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Figure 1: Images of (|47|) . with a = 0.5, for lines of constant radius and constant 
azimuth on the unit disk. 



5.1 The planar case 



For our variables, we replace a point x £ Bd with (x, y), and we replace a point 
s6f! with (s, t). Define the mapping $ : B — > f2 by (s, i) = $ (x, y), 



s = a; — y + ax 
t = x + y 



(47) 



with < a < 1. It can be shown that $ is a 1-1 mapping from the unit disk B. 
In particular, the inverse mapping W : £1 — > _B is given by 



1 r 

a 

1 r 



1 + VT+~a (s + i) 
at - (-1 + i/l + a(s-M)) 



(48) 



In Figure [H we give the images in of the circles r = j/10, j = 1, . . . , 10 and 
the azimuthal lines 9 = jir/10, j — 1, . . . , 20. 

The following information is needed when implementing the transformation 
from —Ait + "fu = / on CI to a new equation on £?: 



D$ = J(aj,y) = 
det (J) = 



l + 2ax -1 
1 1 

2(1 + ax) 
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Figure 2: The function u(s, t) of f50]) 



J(x)' 



1 1 



2(1 + ax) V - 1 l + 2ax 



j t / \-i i/ \-t 1/1 ax 

A = J [Xj J (X) — 7T _ 2 2 , r, 

w w 2(1 + ax) V a:E 2a x + ^ ax + 

The latter are the coefficients needed to define A in (fl"4"|) . 
We give numerical results for solving the equation 

- &u{s 1 t)+e s - t u{s,t) =f(s,t), (s,t)efi (49) 

As a test case, we choose 

u(s,t) = e - " 2 cos (7rt) , (s,t)£fl (50) 

The solution is pictured in Figure [2j To find f(s,t), we use (|49l) and ([50)) . We 
use the domain parameter a — 0.5, with 51 pictured in Figure [TJ 

Numerical results are given in Table[T]for even values of n. The integrations 
in (|37[1 were performed with (|42[) : and the integration parameter q ranged from 
10 to 30. We give the condition numbers of the linear system (f3"Tl) as produced 
in Matlab. To calculate the error, we evaluate the numerical solution and the 
error on the grid 



$ (ZijtUij) = $ i r i COS 6j , 7*i SU10J-) 

10' 10 



( r i> 0j) = ( 7^> 777 ) i i = 0,l,...10; .? = 1,...20 
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Table 1: Maximum errors in Galerkin solution u. 



n 


N n 


UnW^ 


cond 


n 


N n 


[|w-^IL 


cond 


2 


6 


9.71£- 1 


14.5 


14 


120 


3.9QE - 5 


6227 


4 


15 


2.87.E-1 


86.1 


16 


153 


6.37S-6 


10250 


6 


28 


5.85E- 2 


309 


18 


190 


8.20£- 7 


15960 


8 


45 


1.16-B- 2 


824 


20 


231 


9.44.B - 8 


23770 


10 


66 


2.26.E-3 


1819 


22 


276 


1.06-E-8 


34170 


12 


91 


2.8LE-4 


3527 


24 


325 


1.24E-9 


47650 




Figure 3: Errors from Table Q] 



The results are shown graphically in Figure [31 The use of a semi-log scale 
demonstrates the exponential convergence of the method as the degree increases. 

To examine experimentally the behaviour of the condition numbers for the 
linear system (|37p. we have graphed the condition numbers from Table Q] in Fig- 
ured] Note that we are graphing N% vs. the condition number of the associated 
linear system. The graph seems to indicate that the condition number of the 
system (|3"T)l is directly proportional to the square of the order of the system, 
with the order given in (|38[). 

For the Poisson equation 

-Au(s,t) = f(s,t), ( s ,t)efl 

with the same true solution as in (|50p , we use the numerical method given in 53 
The numerical results are comparable. For example, with n — 20, we obtain 
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Figure 4: Condition numbers from Table [T] 



||tt — Unlloo = 9.90 x 10 8 and the condition number is approximately 14980. 
5.2 The three dimensional case 

To illustrate that the proposed spectral method converges rapidly, we first use 
a simple test example. We choose the linear transformation 

(X\ - 3x 2 \ 
2xi + x 2 
Xl + x 2 + x 3 J 

so that -Bi(O) is transformed to an ellipsoid f2i; see figure [3] For this transfor- 
mation D<&\ and J\ = det(D$i) are constant functions. For a test solution, we 
use the function 

u(s) = sie S2 sin(s 3 ) (51) 

which is analytic in each variable. 

Table [5] shows the errors and the development of the condition numbers for 
the solution of ([T|) on Oi . The associated graphs for the errors and condition 
numbers are shown in figures [5] and [71 respectively. The graph of the error is 
consistent with exponential convergence; and the condition number seems to 
have a growth proportional to the square of the number of degrees of freedom 
N n . 

Next we study domains which are star shaped with respect to the origin, 
n 2 = {x G M 3 | x = r(p,6,(f>), < p < R(6, 4>)}. (52) 
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Figure 5: The boundary of f2i 



See (|45|) for the definition of T, and R : § 2 — > (0,oo) is assumed to be a 
C°° function. In this case we can construct arbitrarily smooth and invertible 
mappings <f> : -Bi(O) — > f^2 as we will show now. First we define a function 
* : [0, 1] -> [0, 1] 

o, o< P <i 



[r-iP-hr; Hpi" (53) 

the parameter e s G N determines the smoothness of f E C e » [0, 1]. For the 
following we will assume that R(9, </>) > 1, for all 9 and 0; this follows after an 
appropriate scaling of the problem. With the help of t we define the function 
R which is monotone increasing from to R(9, </>) on [0, 1] and equal to the 



Table 2: Maximum errors in Galerkin solution u. 
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Figure 6: Errors from Table [2] 



identity on [0, 1/2], 

R(p,6,4>) :=t(p)R(9,<f > ) + (l-t(p))p 

Because 

£-R(p, 0, 0) = t'(p)(R{6, 4>) - P) + (1 - t( P )) > 0, p e [0, 1] 
dp 

the function S is an invertible function of p of class C 63 " 1 . The transformation 
$ 2 : Si (0) — * f2 2 is defined by 

$ 2 (x) := T(R(p, 9, 0), 0, 0), x = T(p, 0, 0) G Si(0) 

The properties of R imply that <E> 2 is equal to the identity on Si(0) and the 
outside shell 2?i(0) \ Si (0) is deformed by <J> 2 to cover Q 2 \ Si (0). 
For a test surface, we use 

R(0, 4°) = 2 + I cos(20) sin(6>) 2 (7 cos(6>) 2 - 1) (54) 
e s = 5; 

see figures [8]|9] for pictures of <9f2 2 - For our test example, we use u from (|5l"j) . 

The term cos(2</>) sin(6') 2 (7cos(6') 2 — 1) is a spherical harmonic function which 
shows R 6 C°°(§ 2 ), and the factor 3/4 is used to guarantee R > 1. For the 
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Figure 7: Conditions numbers from Table [2] 



transformation <j> 2 we get $2 € C 4 (-Bi(0)), so we expect a convergence of order 
0(n~ 4 ). Our spectral method will now approximate u o <£> 2 on the unit ball, 
which varies much more than the function in our first example. 

We also note that one might ask why we do not further increase e s (see ((55)) ) 
to get a better order of convergence. It is possible to do this, but the price one 
pays is in larger derivatives of mo <£> 2 , and this may result in larger errors for the 
range of n values where we actually calculate the approximation. The search 
for an optimal e s is a problem on its own, but it also depends on the solution 
u. So we have chosen e s = 5 in order to demonstrate our method, showing that 
the qualitative behaviour of the error is the same as in our earlier examples. 

The results of our calculations are given in tableau and the associated graphs 
of the errors and condition numbers are shown in figures [TOl and [TTl respectively. 
The graph in Figure [TT] shows that the condition numbers of the systems grow 
more slowly than in our first example, but again the condition numbers appear 
to be proportional to N% . The graph of the error in Figure [TU] again resembles 
a line and this implies exponential convergence; but the line has a much smaller 
slope than in the first example so that the error is only reduced to about 0.02 
when we use degree 16. What we expect is a convergence of order 0(n~ 4 ), but 
the graph does not reveal this behavior in the range of n values we have used. 
Rather, the convergence appears to be exponential. In the future we plan on 
repeating this numerical example with an improved extension $ of the boundary 
given in (f5"i|) . 

When given a mapping ip : dB — ► dfl, it is often nontrivial to find an 
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Figure 8: A view of dfl2 



extension <E> : B — > f2 with $|q B = >p and with other needed properties. For 

onto 

example, consider a star-like region f2 whose boundary surface dQ is given by 

P = R{0,4>) 

with R : S 2 — ► c?f2. It might seem natural to use 

&(p,Q,<f>) =pR(d,(j>), < p < 1, 0<6»<tt, O<0<2tt 

However, such a function $ is not continuously differentiable at p = 0. We 
are exploring this general problem, looking at ways of producing $ with the 
properties that are needed for implementing our spectral method. 
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